General analysis of mathematical models for bone remodeling 
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Bone remodeling is regulated by pathways controlling the interplay of osteoblasts and osteoclasts. 
In this work, we apply the method of generalized modelling to systematically analyse a large class of 
models of bone remodeling. Our analysis shows that osteoblast precursors can play an important role 
in the regulation of bone remodeling. Further, we find that the parameter regime most likely realized 
in nature lies close to bifurcation lines, marking qualitative changes in the dynamics. Although 
proximity to a bifurcation facilitates adaptive responses to changing external conditions, it entails 
the danger of losing dynamical stability. Some evidence implicates such dynamical transitions as a 
potential mechanism leading to forms of Paget's disease. 
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I. INTRODUCTION 

Bone is a complex tissue that is being repaired and 
rebuilt continuously throughout an individual's life. 
The process of bone remodeling consists of two sub- 
processes: the resorption of old bone and the forma- 
tion of new bone. In the past decades it has become 
clear that, at the cellular level, bone remodeling de- 
pends on the interplay between two different types 
of cells, osteoclasts and osteoblasts. The former are 
multinuclear cells of hematopoietic origin that resorb 
bone, while the latter are mononuclear cells of mes- 
enchymal origin that fill the gaps left by osteoclasts 
with newly formed bone tissue [lj [2] . 

From a theoretical perspective, bone remodeling is 
interesting because biological requirements seem to 
pose contradictory demands. On the one hand, the 
system must show robustness with respect to natu- 
rally occurring fluctuations. On the other hand, the 
system must show adaptivity to relevant changes in 
the external conditions that require an increased or 
decreased rate of bone remodeling. 

In humans, the process of bone remodeling is reg- 
ulated by several autocrine and paracrine factors to 
maintain the balance of bone. In particular, it has 
been discovered that a signaling pathway involving 
the Receptor Activator of NF-kB (RANK), its ligand 
RANKL, and the cytokine osteoprotegerin (OPG) 
play an important role in the regulation of bone re- 
modeling [Sill]. For osteoclasts to mature, it is neces- 
sary that RANKL, expressed by cells of osteoblastic 
lineage, attaches to RANK, expressed on cells of os- 
teoclastic lineage. This process is regulated by the de- 
coy receptor OPG, which is expressed by osteoblastic 
cells and inhibits the differentiation of osteoclasts by 
binding to RANKL and thus sequestering it. Another 
important regulator, the cytokine TGF/3, is known to 
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influence both osteoclasts and osteoblasts Over- 
or underexpression of TGF/3 and the protagonists of 
the RANKL pathway is related to several diseases of 
bone, such as osteoporosis and Paget's disease of bone 

[MS]- 

At a particular site, osteoblasts and osteoclasts 
move through bone in a group, remodeling tissue on 
its way. Such a collection of cells is known as a basic 
multicellular unit (BMU) [10] , The dynamics of a sin- 
gle BMU is difficult to model because spatial aspects 
become crucial. In this study, mathematical models 
of ordinary differential equations (ODEs) are stud- 
ied, which is an approximation that has often been 
applied in earlier studies. The justification for the 
omission of spatial effects is that an average is taken 
over many BMUs, thus analysing systemic properties 
of bone remodeling. 

Mathematical models of bone remodeling were 
studied in various earlier works. In particular, a mini- 
mal model consisting of two ODEs was constructed in 
Ref. [TT] . In this model, the terms describing the reg- 
ulation were assumed to be power laws. Thereby, the 
effects of paracrine and autocrine factors were con- 
densed into power law exponents. A more detailed 
model was formulated in [T^HLi] which incorporated 
three dynamic variables and made use of Michaelis- 
Menten kinetics. 

The earlier works pointed out that under physi- 
ological conditions and in absence of external stim- 
uli, the system should reside in a steady state, where 
the numbers of osteoclasts and osteoblasts remain ap- 
proximately constant in time. For the system to re- 
main close to the steady state, the state has to be 
dynamically stable, so that the system is driven back 
to the steady state after small perturbations. At the 
same time it is desirable that the stationary densi- 
ties of osteoblasts and osteoclasts react sensitively 
to external influences, communicated through the 
signalling molecules. Mathematically speaking this 
means that the system should be robust against fluc- 
tuations of the variables, but sensitive to changes in 
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the parameters. In dynamical systems, the strongest 
response of steady states to paramater change is of- 
ten found close to bifurcations - critical thresholds 
at which the stability to perturbations is lost. There- 
fore, it is intuitive that there should be some tradeoff 
between dynamical stability and responsiveness. It is 
thus possible that the physiological state of the bone 
remodeling system is characterized by parameter val- 
ues close to a bifurcation point. 

An intriguing possibility, raised in [11] , is that some 
diseases of bone might have their cause not in a shift 
of the steady-state concentrations, but in a bifurca- 
tion, in which the stability of the steady state is lost. 
Dynamical systems theory has established a large va- 
riety of powerful tools for detecting and analyzing 
bifurcations. If a given disease were found to be re- 
lated to a bifurcation phenomenon, this arsenal of 
tools could be utilized for understanding the causes 
and consequences of the disease. 

In dynamical systems, the dynamics close to steady 
states are governed by the Jacobian matrix [15l [16] . 
In a Af-dimensional (i.e., N- variable) system the Ja- 
cobian has N eigenvalues which are either real, or 
form complex-conjugate eigenvalue pairs. Local bi- 
furcations occur when a change of parameters causes 
one or more eigenvalues of the Jacobian matrix to 
cross the imaginary axis, so that an eigenvalue with 
a negative real part becomes an eigenvalue with a 
positive real part. This can occur in two fundamen- 
tal scenarios: In the case of a saddle-node bifurcation, 
a real eigenvalue crosses the imaginary axis and be- 
comes positive. This bifurcation typically occurs at 
a threshold at which steady states collide and van- 
ish, which can lead to abrupt transitions in the sys- 
tem. In the case of the Hopf bifurcation, a complex 
conjugate pair of eigenvalues crosses the imaginary 
axis, acquiring positive real parts, while the stabil- 
ity of the steady state is lost. We can distinguish 
between two types of Hopf bifurcations: In a super- 
critical Hopf bifurcation, a stable limit cycle is born, 
leading to small-amplitude sustained oscillations. In 
a subcritical Hopf bifurcation, an unstable limit cycle 
that coexists with a stable steady state vanishes, typ- 
ically leading to a catastrophic loss of stability and, 
often, large-amplitude oscillations |16) . 

In the present work, we explore a large class of 
mathematical models for the regulation of bone re- 
modeling with respect to their bifurcation proper- 
ties. To this end, we apply the method of general- 
ized modeling [T?H2"D] which allows analyzing models 
in which the reaction kinetics is not restricted to spe- 
cific mathematical functions. Thereby, generalized 
models can provide a broad overview of the dynamics 
of the system, which facilitates the choice of specific 
models. More importantly, the analysis of the gener- 
alized model reveals dynamical instabilities that are 
potentially related to pathologies of bone remodeling. 

The models proposed in the present paper general- 



ize and extend findings from earlier specific models. 
Specifically, we find that among the previously dis- 
cussed scenarios for regulatory interactions the one 
that is known to yield the highest responsiveness 
leads to steady states close to bifurcations. Cross- 
ing these bifurcations causes trajectories to leave the 
dynamical regime of healthy bone remodeling and can 
possibly be related to diseases of bone. We show that 
in two- variable models with the structure proposed in 
pTj , stability of the steady state requires that the ac- 
tions of OPG dominate over those of RANKL, while 
three-dimensional models [TH H3] can also be stable 
without that assumption. These results suggest that 
osteoblast precursors have an important impact on 
the dynamics and should be taken into account ex- 
plicitly in models. 



II. MODEL CONSTRUCTION 

A mathematical model capable of describing the 
process of bone remodeling must take into account 
the concentrations of active osteoblasts, B, and of 
active osteoclasts, C. We start by discussing a mini- 



mal model of these two dynamic variables in Sec. II A 



Because the minimal model potentially oversimplifies 
the problem by ignoring the populations of precur- 
sor cells, especially in the case of osteoblasts, a more 
detailed model that accounts for osteoblast behav- 
ior at different stages of maturation is introduced in 
Sec. HTBl 



A. Two-variable model 

In the minimal two-variable model, we assign to 
both state variables a gain term (F,H, respectively) 
describing the maturation of new cells from a pool of 
precursors, and a loss term (G,K) describing the re- 
moval of cells due to death or further differentiation. 
This leads to the basic equation system 



dt 
d 

dt 



B = F(B,C)-G(B,C) 
C = H{B,C)- K{B,C) 



(1) 



In the following, we assume that the functions 
F(B,C),G(B,C),H(B,C) and K{B,C) are positive 
and continuously diffcrentiable, but do not restrict 
them to specific functional forms. 

Our aim is to derive mathematical conditions on 
the stability of all possible steady states in all mod- 
els of the form introduced in ([I]). As we show below, 
the threshold parameter values at which the stability 
changes (i.e., the bifurcation points) can be expressed 
as a function of parameters, which have a clear biolog- 
ical interpretation in the context of bone remodeling. 
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In any particular steady state we can denote the 
number of osteoblasts and osteoclasts as B* and C* , 
respectively. We can then define normalized variables 

Similarly, we define a set of normalized functions 

F(B,C) G(2?,C) 
/(M = ip* i p* 9(b,c) = 



p*(B*,C*)- 



G*(B*,C*Y 



H(B,C) K(B,C) 
h(b > c) = H*(B*,C*y k( > b > c)= K*{B*,C* 



(3) 



Using these definitions, the system can be written as 
d b = ai (/(6,c) - g(b,c)) 

(4) 



dt 
d 



dt 



c = «2 (M^j c) — k(b, c)) . 



where 



and 



F*(B*,C*) G*(B* : C*) 

a l = = ^ ( 5 ) 



B* 



H*(B*,C*) K*(B* 7 C*) 

"2 = 7^2 = ■ I 6 ) 



C* 



c* 



The second equalities in Eq.([5j| and Eq.([6]) hold be- 
cause gain and loss terms have to balance in a steady 
state. 

In the new coordinates, the formerly unknown 
steady state is located at (b, c) = (1,1). As men- 
tioned in the introduction, this steady state is said to 
be asymptotically stable if the system returns to it 
after a sufficiently small perturbation [HI [16] . This 
is the case if all eigenvalues of the Jacobian matrix 
have negative real parts. In the present model the 
Jacobian can be written as 



"i \ f f h - g h f c - g c 
a 2 I \hb - kb h c - k c 



(J) 



Here we used Roman subscripts to indicate partial 
derivatives. For instance, /b is defined as 



fb 



' db 6=l,c 

d (InF) 



B* dF 
~F*dB 



B=B',C=C* 



d (InS) 



(8) 



B=B',C=C* 



So far we succeeded in constructing the Jacobian 
matrices corresponding to steady states in a very gen- 
eral class of models. We emphasize that we did not 
have to assume that there was only one steady state 
in the system. In the general case where more steady 
states exist, the formal derivation of the Jacobian ap- 
plies to all steady states in all models within the class 



considered here, whereas the quantities appearing in 
the Jacobian matrix generally differ between steady 
states. Although these quantities, such as /b, are in 
general unknown, they do not depend on the dynami- 
cal variables and can therefore be treated as unknown 
parameters with the same right as the parameters 
that are introduced in conventional models. Just as 
conventional parameters, the generalized parameters 
appearing in the Jacobian have a well-defined inter- 
pretation in the context of the model. In the remain- 
der of this subsection we discuss the interpretation in 
detail. 

The parameters ot\ and a 2 are defined as ratios be- 
tween a flux and a concentration and thus have the 
dimension of an inverse time. They represent the re- 
spective time scales of the two coupled differential 
equations and can be interpreted as the inverse life- 
time of the respective cell type. Since the average life 
span of osteoblasts (~3 months) exceeds the life span 
of osteoclasts («2 weeks) by a factor close to 6 [T] , it 
is reasonable to assume that a\ja 2 « 1/6. Since the 
scale by which time is measured is arbitrary, we are 
free to fix a.\ = 1. 

As can be seen from Eq. ([8]), the remaining param- 
eters in the Jacobian are logarithmic derivatives of 
the original functions. We denote these parameters 
as elasticities, using a term from Metabolic Control 
Analysis [21] . In the simple case of a linear functional 
dependency, the corresponding elasticity is exactly 1. 
In the case of a power-law function depending on a 
variable x, f(x) — ax b , the elasticity is the exponent 
b. We note that in the model proposed in Ref. [H] 
all processes were modeled as power-laws. There- 
fore the Jacobian that is derived in this earlier paper 
is mathematically identical to Eq. ([7|. Nevertheless 
the Jacobian derived in the present work describes a 
larger class of models, in which processes are mod- 
eled by arbitrary positive functions. For modeling 
dependencies that saturate for high concentrations, a 
Michaelis-Menten-type function is often assumed. In 
this case, the corresponding elasticity is confined to 
the interval [0, 1]. It is close to 1 in the initial region 
of steepest slope and approaches in the regime close 
to saturation. 

If the exact functional form is not known, it is pos- 
sible to assign a range of plausible values to the elas- 
ticities that is based on biological knowledge. Elas- 
ticities that are associated with an activating influ- 
ence (positive feedback) are positive, whereas elas- 
ticities associated with an inhibiting influence (neg- 
ative feedback) are negative. For instance, many 
reasonable feedback functions, such as the Michaelis- 
Menten function from enzyme kinetics, first grow lin- 
early with the argument but then approach satura- 
tion for large values of the argument. For such a 
function, the corresponding elasticity parameter is re- 
stricted to the interval [0, 1] |22) . 

In the following we assume that the osteoblasts' 
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lifetime is not affected by additional regulators. 
Therefore, the decay term of osteoblasts is linear in b, 
and is independent of c. This translates into g c — 
and gb = 1. Likewise, the decay term of osteo- 
clasts is not influenced by osteoblasts, correspond- 
ing to fcb = 0. Moreover, there is no evidence for 
strongly nonlinear autocrine regulation of osteoblasts 
(fb ~ 0). Nevertheless, we analyze the bifurcations 
with respect to fb in order to determine if a posi- 
tive or negative feedback mechanism would affect the 
dynamics. 

The parameters f c , h c and k c depend on the growth 
factor TGF/3, an important regulator in bone remod- 
eling [5]. When osteoclasts resorb bone tissue, TGF/3 
is released into the bone matrix, where it facilitates 
the differentiation of osteoblast progenitors to active 
osteoblasts, leading to f c > 0. 

The autocrine roles of TGF/3, described by the pa- 
rameters h c and k c , are less clear: In vitro experi- 
ments have led to contradictory results on the influ- 
ence of TGF/3 especially on osteoclasts, finding both 
activation and repression [5] . Results depend strongly 
on the experimental setup, such as TGF/3 concentra- 
tion or whether isolated cultures or co-cultures are 
used. According to a current view [5], TGF/3 in- 
directly acts as a repressor by interaction with the 
OPG/RANKL/RANK pathway in co-cultures with 
osteoblasts. In isolated cultures however, TGF/3 ac- 
tivates and sustains osteoclasts. For covering both 
cases, we allow the corresponding parameter, h ci to 
be positive or negative. 

Without any feedback mechanisms, one would as- 
sume that the decay term of osteoclasts, k Cl were 
equal to one. However, the apoptotic decay of osteo- 
clasts has been reported to both be promoted [25H25] 
and suppressed [H] [27J by TGF/3, corresponding to 
k c > 1 and k c < 1, respectively. Based on these 
conflicting experimental results, we assume that the 
corresponding elasticity k c is positive but do not as- 
sume a specific value. 

Because of the form of the Jacobian, Eq. only 
the difference m c = h c — fc c is important for the sta- 
bility analysis. Therefore, the effects of autocrine reg- 
ulation in the production and decay terms of osteo- 
clasts can be covered by a single parameter. 

Finally, the effects of the RANKL /RANK/OPG 
system are condensed into the parameter hb- This 
parameter can be negative or positive, depending on 
whether the repressive influence of OPG or the acti- 
vation of RANK dominates. 

In summary, the considerations above lead us to 
the Jacobian 
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FIG. 1: Schematic sketch of the two-variable model. 
Osteoblasts (OB) influence osteoclasts (OC) via the 
RANKL /RANK/OPG pathway, while the TGF-/3 path- 
way exerts a positive feedback from osteoclasts to both 
osteoclasts and osteoblasts. 



B. Three-variable model 



It can be argued that the two variable model pro- 
posed above is oversimplified because the dynamics 
of precursor populations is neglected. In particular, 
one might miss important information by describ- 
ing the population of osteoblasts with a single dy- 
namic variable when osteoblasts have different prop- 
erties at different stages of maturation [351 US]- A 
model with a different structure has been proposed 
in [T2] and was subsequently extended in [T3 HI] . 
In this model, cells of osteoblastic lineage are rep- 
resented by two dynamic variables, responding os- 
teoblasts (ROBs), R and active osteoblasts (AOBs), 
B. ROBs are committed to the osteoblastic lineage 
and interact with osteoclasts but are not yet func- 
tional osteoblasts. There are two important reasons 
for distinguishing AOBs and ROBs: First, there is 
experimental evidence that osteoblastic cells express 
RANKL and OPG differently at different stages of 
maturation, where at later stages the ratio of RANKL 
to OPG seems to decrease [2E1I2S]- Second, TGF/3, 
which is released and activated by osteoclasts, acti- 
vates osteoblast differentiation only at early stages of 
differentiation, whereas it seems to enlarge the pool 
of responding osteoblasts by inhibiting further differ- 
entiation into active osteoblasts [5]. 



with three remaining free parameters. 



The structure of the three-dimensional model, 
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FIG. 2: Schematic overview of the three-variable model, 
in which the dynamics of responding osteoblasts (ROB), 
active osteoblasts (AOB) and osteoclasts (OC) is de- 
scribed. The feedback mechanisms, mediated by the 
RANK/RANKL/OPG-pathway and by TGF/3, are in- 
scribed in the diagram in the form of arcs with arrows. 
The straight arrow from ROB to AOB indicates a flow of 
biomass due to differentiation of ROBs. 



shown in Fig. [2] translates to the equations 
(1 



dt 
d_ 

dt 
d_ 

dt 



R = S{C)-T(R,C) 
B = T(R,C)-U(B) 
C = V(B,R)-W(C), 



where the two terms in each line again correspond 
again to gains and losses of the population of the 
respective cell type. The functional dependence in 
these equations is motivated by the biological pro- 
cesses that are included in the model. We explain 
these processes in more detail after formally con- 
structing the Jacobian. 

Performing the normalization procedure that was 
outlined in the previous section and in [T7], we can 
describe the structure of the model with the normal- 
ized equations 



d 

dT 

dt 
d 

dt' 



an (s(c) - t(r, c)) 
a 2 {t(r, c) - u(b)) 
a 3 (v(b,r) - w(c)) 



where, in analogy to our treatment of the two- variable 
model, the lower-case variables and functions denote 
the normalized quantities and ot\, 0:3 are the char- 
acteristic timescales of ROB, AOB, and osteoclast 
turnover. 

In analogy to Eq. Q, the Jacobian for the three- 
variable model can now be written as 



(10) 



A summary of the elasticities occurring in the Jaco- 
bian and the ranges we assign to them is given in 
Table 1. 
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The elasticities s c , t c and w c describe the nonlin- 
earities in the TGF/3 pathway. This pathway stabi- 
lizes the reservoir of AOBs both by promoting the dif- 
ferentiation of osteoblast progenitors to ROBs, lead- 
ing to s c > and by inhibiting the further differen- 
tiation of ROBs to AOBs, leading to t c < 0. In our 
model, we restrict s c to the interval [0, 1] and t c to 
[—1,0]. This range includes the choice of Hill func- 
tions with exponents equal to 1 that were used in 
earlier models such as [14]. The nature of autocrine 
regulation of osteoclasts has not been ultimately clar- 
ified. Therefore we restrict w c to [0.5,1.5]. This 
range is centered around w c = 1, because without 
any additional feedback, a linear decay term would 
be expected. In particular, the parameter is smaller 
(greater) than one if the additional feedback is neg- 
ative (positive). We note that the functional forms 
that were assumed in Ref. 12, 13] lead to superlinear 
decay (w c > 1). 

The regulation of osteoclasts by cells of osteoblas- 
tic lineage is mediated by the RANKL/RANK/OPG 
pathway. Depending on the ratio between RANKL 
and its decoy receptor OPG, the corresponding elas- 
ticities v T and Vb can be either positive or negative. 
This includes all possible combinations of RANKL 
and OPG expression at responding osteoblasts or ac- 
tive osteoblasts. In particular, two fundamentally dif- 
ferent scenarios are described in the literature, which 
are for instance discussed as models Ml and M2 in 
Ref. [13] • These scenarios are characterized in the 
general model by 

1. v T < 0, Vb > 0. OPG is expressed by respond- 
ing osteoblasts, RANKL is expressed by active 
osteoblasts. 

2. v T > 0, Vb < 0. RANKL is expressed by re- 
sponding osteoblasts, OPG is expressed by ac- 
tive osteoblasts. 

Intermediate situations with a differential expression 
of OPG and RANKL without the assumption of ex- 
clusive expression are also covered by our description 
(e.g. v r > v h > 0). 



III. RESULTS 

In this section, we show results for the two models 
that were introduced in the preceding section. The 
bifurcations in both models can be computed analyt- 
ically from the Jacobians. In the two-variable model, 
the results of the bifurcation analysis depending on 
three parameters can be directly visualized and un- 
derstood intuitively. However, in the three-variable 
model the larger number of parameters complicates 
gaining intuitive understanding. For our initial explo- 
ration of the generalized model we therefore resort to 
a numerical procedure which provides results that are 
more easily interpretable. 
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A. Bifurcation analysis of the two-variable 
model 

In any system, a necessary condition for a saddle- 
node bifurcation is det J — 0, guaranteeing a zero 
eigenvalue. For the Jacobian derived in Eq. |7| it 
follows that 



RANKL dominates over repression by OPG, desta- 
bilizes the steady state in a saddle-node bifurcation 
(lower region of Fig. [3| , while a negative hy, has a sta- 
bilizing effect. In the model, stability is therefore pro- 
moted when OPG dominates over RANKL and the 
effective action of the RANKL pathway is inhibiting. 
It is not clear whether this requirement is fulfilled in 



m c (f h - 1) - f c h h = 



(11) 



has to be satisfied at a saddle-node bifurcation (m c = 
he &c) ■ 

In the Hopf bifurcation, there are two conjugate 
eigenvalues with zero real part. In a two-dimensional 
system, this means that the eigenvalues add up to 
zero and the trace of the Jacobian vanishes (Tr J = 0) 
so that 



-(/b-l) + 



0. 



Qf2 



(12) 



is a necessary condition for the Hopf bifurcation. Ad- 
ditionally, the inequality det J > must be fulfilled. 




FIG. 3: Bifurcation diagram for the 2- variable model, de- 
pending on /b, m c and —f c hb. Each combination of the 
parameters in the three-dimensional volume corresponds 
to the steady state in a particular model. There are two 
distinct bifurcation surfaces that divide the regions where 
steady states are stable (SS) from regions where steady 
state unstable (US). The red surface is formed by Hopf 
bifurcation points, whereas the blue surface is formed by 
saddle-node bifurcation points. 

Because the ratio of timescales cci/a.2 is known, the 
bifurcation conditions, Eq. ( 11 ) and Eq. ( 12 ), depend 



only on different combinations of the three parame- 
ters, /b, m c and —f c hb (for stability, only the product 
appearing in the Jacobian is important). The stabil- 
ity of all steady states in the whole class of models can 
therefore be visualized in the single three-parameter 
bifurcation diagram, displayed in Fig. [3] 

The figure shows that parameter regimes of insta- 
bility can be reached both via a Hopf or a saddle-node 
bifurcation. Because f c > 0, we see that a high value 
of /ib, corresponding to the case where activation by 



The parameter regime close to the Hopf bifurca- 
tion can only be reached when m c « 0, i.e., when 
the autocrine feedback of the osteoclasts acts in an 
activating way as a countereffect to the linear contri- 
bution of the decay term. 

We further note that the bifurcation diagram dif- 
fers from Fig. 4 a) in Ref. [TT], in which the saddle- 
node bifurcation surface seems to be independent 
from /b (called 322 there), which is incompatible with 
the form of Eq. [TT] 

In the generalized model we cannot determine 
wether the Hopf bifurcation is subcritical or super- 
critical without making further assumptions. We 
therefore consider the model proposed in [TTJ as a 
specific example of the more general class considered 
here. In our notation this model can be written as 



±B = A b Bf*Cf° - D b B 
4(7 = A c B h »C h < - D C C, 



(13) 
(14) 



where Ab,D},,A c and D c are rate constants The Hopf 
bifurcation condition in this model is 



(h c - 1) + (/b - 1) = 0. 



(15) 



This equation is equivalent to Eq. A5 in [TT] . 

However, the very same condition guarantees that 
the flow is Hamiltonian, i.e., that a function of the 
dynamic variables exists which is conserved on all tra- 
jectories. We found this function to be 



H(B, C) = -J^—B^-^C 1 -^ 
76 — 1 



A c . 



h-fb + 1 



B 



h b -f b + l 



A h 



fc - h c 



(16) 



It can easily be verified by d irec t differentiation that 
under the condition of Eq. (15), ^ = 0. The ac- 
tual Hamilton equations are fulfilled after the coordi- 
nate transformation p(B) = jz~rB 1 ~^ b and q(C) — 



C 



l-hc 



It follows that no limit cycle with a defined ampli- 
tude is born in the Hopf bifurcation and the steady 
state is a center. The Hopf bifurcation is neither 
subcritical nor supercritical, but is just at the brink 
between the two alternatives. This structural insta- 
bility is caused by a symmetry in the model. Sus- 
tained oscillations occur only exactly at the bifurca- 
tion point and the amplitude of the oscillations will 
depend strongly on initial conditions. 
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The symmetry that causes the degenerate behavior 
of the model is broken if the autocrine regulation of 
the osteoclast removal is taken into account (fc c 7^ 1), 
as it has been assumed in other models jT^j- Both 
for negative and positive autocrine feedback, the sys- 
tem is no longer Hamiltonian at the Hopf bifurca- 
tion, and depending on the actual parameters the bi- 
furcation is subcritical or supercritical. We verified 
numerically that close to supercritical Hopf bifurca- 
tions, stable limit cycles with a well-defined ampli- 
tude can exist and sustained oscillations are possible 
over a wider range of parameters. These findings im- 
ply that the two-dimensional mathematical model is 
sensitive with respect to the existence of feedback in 
the removal term of osteoclasts.. 

In summary, the bifurcation analysis shows that in 
the two-dimensional model structure, it is necessary 
for stability that the feedback exerted by osteoblasts 
on osteoclasts is effectively inhibiting, which is the 
case when the repressing effects of OPG dominate 
over the activating effects of RANKL. Otherwise, 
the steady state under consideration is unstable and 
therefore cannot correspond to a physiological equi- 
librium. Furthermore, we showed that a Hopf bifur- 
cation exists close to the realistic parameter regime. 
A pathological shift in the parameters may therefore 
drive the system over the Hopf bifurcation, which will 
typically lead to stable sustained oscillations of osteo- 
clast and osteoblast numbers if the osteoclast removal 
rate increases at least weakly with the number of os- 
teoclasts. 



B. Bifurcation analysis of the three-variable 
model 
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FIG. 4: Effect of parameters on local dynamics. Using an 
ensemble of 10 6 randomly drawn steady states, the his- 
tograms show the percentage of randomly drawn states 
that are unstable (red crosses) and the fraction of ran- 
domly drawn states which are unstable and have lead- 
ing complex eigenvalues, indicating oscillatory instabili- 
ties (green circles). Each panel shows the effect of one of 
the elasticities, while averaging out the effect of the other 
parameters. An ascending curve signifies that for values 
at the top of the range of the respective parameter, steady 
states are more likely to be unstable than for low values 
of the parameter, while a descending curve signifies the 
opposite. 



In the three variable model we cannot visualize all 
factors impinging on stability in a single diagram. We 
therefore use a random sampling procedure |18j for 
gaining a first impression of the effect of the various 
parameters. This analysis is then combined with a 
bifurcation analysis of three-dimensional subspaces of 
the larger parameter space. 

For the three-variable model, we performed a sta- 
tistical sampling sampling analysis by creating N — 
10 7 random parameter sets. In each set, we assigned 
to all parameters random values that were drawn 
from uniform distributions, using the intervals de- 
fined in Table H 

We then determined the stability of the steady 
state in each sample by numerically computing the 
eigenvalues of the Jacobian. Based on this ensemble 
of random steady states, we then statistically analyse 
the relations between the parameters and stability by 
creating histograms for each parameter, showing the 
percentage of unstable steady states in the whole en- 
semble as a function of the parameter value. The 
results of this analysis are shown in the histograms 



Parameter 


Interpretation 


Range 


Sc 


activation of ROB production 


[0,1] 


tr 


ROB decay, linear in r 


1 


tc 


repression of ROB decay 


[-1,0] 


Ub 


AOB decay, linear in b 


1 


Vi 


action of ROBs on OC 


[-1,1] 


Vb 


action of AOBs on OC 


[-1,1] 


W c 


activation of OC decay 


[0.5, 1.5] 


TABLE I: 


Parameters in the three- variable model 



in Fig. [4] 

Panel A and B of Fig. [3] show that a strong 
paracrine effect of TGF/3 on osteoclasts (s c >• 
and t c <C 0) destabilizes steady states (i.e., on the 
left side of Fig. [4j^V where s c is small, the propor- 
tion of randomly generated states that are unstable 
is close to zero). The parameter w c has the oppo- 
site effect (Fig. [4]G) : Strong positive feedback of os- 
teoclasts on osteoclast removal stabilizes the steady 
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state. It follows that the paracrine effects of TGF f3 
on osteoblasts that are described by s c and t c destabi- 
lize the steady state, whereas the apoptosis-inducing 
autocrine effects of TGF/3, described by w c , stabilize 
it. 

For v r < 0, we detected few unstable states 
(Fig. |4fD), showing that models in which only OPG 
is preferentially expressed on ROBs usually operate 
from a stable steady state that cannot be destabi- 
lized easily. The other parameter that is related 
to RANKL signaling, Vh, also acts destabilizingly 
at large positive values (Fig. |4E). As noted above, 
there is experimental evidence that OPG is expressed 
stronger on active osteoblasts, while RANKL is ex- 
pressed stronger on ROBs. This implies that the pa- 
rameter regime that is most likely realized in nature 
is characterized by v T > «b, which is also the regime 
in which instabilities occur most frequently. 

In order to investigate the nature of the instabili- 
ties in more detail, we distinguished between unstable 
steady states in which the eigenvalue with the largest 
real part is a real number and those in which it is 
part of a complex conjugate pair. The significance 
of this eigenvalue lies in its effect on the departure 
of the system from the unstable state. Specifically, 
when departing from a state in which the eigenvalue 
with the largest real part has a non-zero imaginary 
part, the system launches into oscillations. 

Figure [4] shows that for most parameters the frac- 
tion of unstable states with a complex leading eigen- 
value changes proportionally to the total fraction of 
unstable states for most parameters. However, very 
different behavior is observed for the parameter v-^ 
(Fig. [4|3): While the fraction of unstable states with 
a real positive eigenvalue increases with increasing 
Wb, the fraction of unstable states with a leading pair 
of complex conjugate eigenvalues decreases with an 
increasing v^. This behavior suggests that the main 
route to instability is via a saddle-node bifurcation, 
but the probability for encountering a Hopf bifurca- 
tion increases with decreasing Vb- 

Now that we have identified the overall impact of 
the parameters on stability, we proceed by investi- 
gating selected parameters in bifurcation diagrams. 
Here, we chose to concentrate on the parameters 
Vh and u r for which the random-sampling analysis 
turned out interesting results, as well as the param- 
eter w c , which captures the different ways in which 
osteoblast removal has been modeled in earlier mod- 
els. 

The bifurcation diagram in Fig. [5] shows that pa- 
rameter sets corresponding to stable steady states are 
characterized by large values of w c and small values of 
v r (upper front of the figure). The section in param- 
eter space that is most likely realized in nature based 
on experimental results is characterized by w c w 1, 
which can be close to a Hopf bifurcation depending 
on the values of t> r and Vh that describe whether OPG 




FIG. 5: Bifurcation diagram of the extended model, de- 
pending on the effects of RANKL/OPG (i>b, v r ) and the 
autocrine effects of osteoblast decay (w c ). The stable pa- 
rameter regime (SS), which is located in the upper front 
part of the diagram, can be lost via a Hopf bifurcation 
(red) or a saddle-node bifurcation (blue). Other parame- 
ters: Ql = 1, Q2 = 1, 03 — 6, s c — 0.8, t c = —0.8. 



and RANK are expressed preferentially on osteoblast 
precursors or active osteoblasts. Moreover, Fig. [5] 
confirms the findings from the random-sampling anal- 
ysis that both large values of v T and small values of 
w c act in a destabilizing way. 

For the parameter v^, the situation is more compli- 
cated: For large values of v\>, the steady state loses its 
stability in a saddle-node bifurcation, whereas stabil- 
ity is lost in a Hopf bifurcation for small values. This 
explains that the stability curve for the parameter 
in Fig. [4}D depends on whether all unstable states or 
only those with an oscillatory instability were taken 
into account. The minimum observed for v r in Fig. [4] 
can be explained by the saddle-node bifurcation sur- 
face replacing the Hopf bifurcation surface as the pri- 
mary source of instability when v r is increased. 

We note that the bifurcation diagram in Fig. [5] 
contains several bifurcations of higher codimension. 
While the bifurcations discussed so far are of codi- 
mension 1, forming two-dimensional planes in a 
three-dimensional bifurcation diagram, bifurcations 
of higher codimension appear as lines or points. The 
Hopf-bifurcation surface ends in a Takens-Bogdanov 
bifurcation of codimension two as it connects to the 
saddle-node bifurcation surface. For low values of 
w c , the Hopf-bifurcation intersects the saddle- node 
bifurcation in a Gavrilov-Guckenheimer bifurcation. 
In the center of the Figure, the Takens-Bogdanov bi- 
furcation and the Gavrilov-Guckenheimer bifurcation 
intersect in a triple point bifurcation. The presence 
of codimension-2 bifurcations can be of relevance for 
applications because they can imply the existence of 
non-local properties such as homoclinic bifurcations 
or chaos. However, a detailed discussion of these im- 
plications is beyond the scope of the present paper. 
Instead, we refer the reader to [I61I30] . 
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FIG. 6: Bifurcation diagram in which all processes con- 
trolled by TGF/3 exhibit the same degree of nonlinear- 
ity. The strength of this nonlinearity is described by the 
parameter c = s c = — t c = tu c — 1. The other bifur- 
cation parameters, w r and v\>, describe the effect of the 
RANKL pathway. In the diagram, the red surface de- 
scribes a Hopf bifurcation, whereas the blue surface de- 
scribes a saddle-node bifurcation. In the front region of 
the diagram, steady states are stable. Other parameters: 
ai = 1, Oii = 1, ct3 = 6 



A different section of the parameter space was 
considered in Ref. [12], where a single Hill function 
with one K m value was chosen for all processes con- 
trolled by TGF/3. In the case of repression, the in- 
verse of this function was used. Translated into the 
framework of generalized modeling, this means that 
s c = — t c = w c — 1 = c and thus a reduction of free 
parameters. The new parameter c changes simultane- 
ously the local nonlinearity of all functions describing 
the TGF/3 pathway. The diagram in Fig.[6]shows that 
an increase of the effective feedback by TGF/3 can 
lead to instability, which was not observed in Fig. [5] 
where only the parameter w c was varied. The bifurca- 
tion properties with respect to the other parameters, 
v r and «b, are consistent with Fig. [5} 



IV. DISCUSSION 

In this paper, we have investigated the dynamics of 
a large class of models for bone remodeling using the 
approach of generalized modeling. Investigating the 
bifurcation behavior of a two-variable model topol- 
ogy, we showed that both saddle-node bifurcations 
and Hopf bifurcations can occur. 

In the two-dimensional model, stability of the 
steady state requires that the effect of OPG dom- 
inates over that of RANKL to make the two- 
dimensional system an adequate model for the pro- 
cess of bone remodeling. We further showed the pos- 
sibility of negative or positive autocrine feedback on 
osteoclasts should be taken into account in models 
because assuming a linear removal rate can lead to 



structurally unstable models. In such models any 
small deviation from the model assumptions can lead 
to qualitatively different behavior. Because the gen- 
eralized model proposed here does not need to assume 
any specific functional form, it avoids such degenera- 
cies that often are caused by an unfortunate choice 
of functional forms in conventional models. 

In the analysis of an alternative model with three 
variables, we combined a random sampling approach 
with a bifurcation analysis of specific subclasses of 
models. The three-dimensional model incorporates 
experimental findings suggesting that RANKL is 
expressed preferentially on responding osteoblasts, 
while its antagonist OPG is mainly expressed on ma- 
tured active osteoblasts. These conditions place the 
system into an area of the bifurcation diagram that is 
close to both saddle-node and Hopf bifurcations. The 
stability analysis therefore shows that in the dynam- 
ical system of bone remodeling, various bifurcations 
not only exist but are located in a parameter space 
supported by experimental findings. 

The main benefit of operating close to a region of 
instability is probably that a stronger adaptive re- 
sponse to external changes of the model parameters 
is possible. Although our modelling approach is not 
designed to study the response of the model to pertur- 
bations directly, it has been shown before by direct 
simulations that conventional models that fall into 
the general class considered here respond strongly to 
perturbations, thus allowing a better functional con- 
trol in bone remodeling |13| . 

Despite the benefits, operating close to a bifurca- 
tion also poses risks to the system. A change in the 
parameters by an external process can shift the sys- 
tem over the bifurcation, so that the stable steady 
state becomes unstable or ceases to exist. It is there- 
fore reasonable to ask wether certain diseases of bone 
can be related to bifurcations, leading to qualitatively 
different dynamical behavior. 

It is known that several diseases of bone are re- 
lated to dysfunctions in the regulation of bone re- 
modeling, among them postmenopausal osteoporosis, 
Paget's disease, osteopetrosis and osteopenia. There 
is some evidence that diseases may lead to qualita- 
tively different dynamical behavior. Periodic activity 
of osteoclasts has been observed in Paget's disease 
of bone [6] and also in vitro [31]. It is conceivable 
that these dynamics are evoked by the transition of 
a steady state to instability in a Hopf bifurcation. 
In the two-dimensional model, Hopf bifurcations are 
most likely caused by increasing the activating au- 
tocrine feedback of osteoclasts. Yet, TGF/3 was found 
to be not related to Paget's disease of bone [31] • In 
the three-dimensional model, however, stability can 
also be lost in Hopf bifurcations by increasing the 
RANKL/OPG ratio, which is in agreement with find- 
ings that the RANKL pathway is involved in Paget's 
disease. In particular, OPG deficiency was reported 
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to be related to juvenile Paget's disease [7]. 

In conclusion, it is still unclear if known diseases 
of bone can be connected to bifurcation phenomena. 
However, the analysis presented here suggests that 
Hopf and saddle-node bifurcations exist close to the 
physiological steady state. We do not claim that spe- 
cific diseases can be related to these bifurcations. Yet, 
a bifurcation occurring in vivo should certainly lead 



to a pathological condition. Therefore it seems very 
plausible that a connection for instance between the 
crossing of a Hopf bifurcation and the onset of Paget's 
disease may exist. If this is indeed confirmed it would 
imply that powerful tools of bifurcation theory and 
related data analysis techniques, can be applied to 
explore the dynamics of the disease. 
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